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Draft version February 5, 2008 

ABSTRACT 

We formulate the problem of the formation of magnetically supercritical cores in magnetically 
subcritical parent molecular clouds, and the subsequent collapse of the cores to high densities, past 
the detachment of ions from magnetic field lines and into the opaque regime. We employ the six-fluid 
MHD equations, accounting for the effects of grains (negative, positive and neutral) including their 
inelastic collisions with other species. We do not assume that the magnetic flux is frozen in any 
of the charged species. We derive a generalized Ohm's law that explicitly distinguishes between flux 
advection (and the associated process of ambipolar diffusion) and Ohmic dissipation, in order to assess 
the contribution of each mechanism to the increase of the mass-to-flux ratio of the central parts of 
a collapsing core and possibly to the resolution of the magnetic flux problem of star formation. The 
results, including a detailed parameter study, are presented in two accompanying papers. 

Subject headings: ISM: clouds - ISM: dust - magnetic fields - MHD - stars: formation - shock waves 
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1. INTRODUCTION: OBSERVATIONAL AND 
THEORETICAL BACKGROUND 

In the Milky Way, stars are observed to form in- 
side molecular clouds, located primarily along spi- 
ral arms. Molecular clouds do not exhibit signifi- 
cant ordered expansion or contraction velocities. They 
are however observed to have su personic but sub- 

Alfve n ic linewidth s. Po larization (iHildebrand et aL\ 

119991 iRao et al. I 119981: iLai et al. I I2002D and Zee- 



man (jCrutcher. Kazes. fc Troland Ill987t ] ICrutcher et al. I 
Il994[ ) observations have revealed that such clouds are 
threaded by ordered magnetic fields, of magnitude suffi- 
cient to support them against their self-gravity. 

In such magnetically supported {magnetically subcrit- 
ical) clouds, magnetically supercritical fragments (or 
cores) form by the process of ambipolar diffusion, which 
allows neutral molecules to diffuse through charged par- 
ticles and magnetic field lines, on a timescale 1-10 Myr 
(depending on the local physical conditions). When self- 
gravity becomes strong enough to overwhelm the mag- 
netic forces, the fragments begin to contract dynami- 
cally. The initial stages of this contraction are isother- 
mal. However, the column density eventually increases 
enough that radiation is trapped, the temperature rises, 
and hydrostatic protostellar cores form in the fragments. 
These late stages of contraction and the formation of hy- 
drostatic protostellar cores are the subject matter of this 
paper. 

The leftover magnetic flux at the end of the early, 
isothermal phase of contraction (n < 10 



10 



cm 



exceeds observed stellar m agnetic fluxes by tw o to 
four orders of magnitude (|Mouschovias I Il987l ) and 
must therefore be redistributed or dissipated dur- 
ing the pre-main-sequence evolution. Approximate 
theoretical estimates concluded that Ohmic dissipa- 
tion is effective in destroying the magnetic flux at 
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densi ties above 10^° cm ^ ([Pneuman fc Mitchell I 
1965'; Nishi. Nakano. fc Umebavashi I Il991h . 

Pneuman & Mitchcff] (|1965l ) considered both am- 
bipolar diffusion and Ohmic dissipation as magnetic 
flux-loss mechanisms and performed detailed calcula- 
tions of the charged species abundances, accounting for 
ionization due to cosmic rays, radioactivity and thermal 
ionization. They estimated the relative importance 
of the flux-reduction mechanisms by comparing their 
timescales to the free-fall timescale, and concluded that 
Ohmic dissipation is more effective than ambipolar dif- 
fusion at densities > 10^° cm~'^. The field is recoupled to 
the matter when, at high enough temperatures, thermal 
ionization raises the degree of ionization sufficiently. 
They estimated the magnetic field of a newly born 

solar-mass star to be ^ 40 Gauss. 

Two decades later iNakano fc Umebavashi I ()198 6a'.*D) 
repeated these calculations using a more refined chem- 
ical model. Through timescale comparisons for a 
series of different geometries (disk- like or spherical), 
they also found that the magnetic field decouples from 
the matter because of Ohmic dissipation. This de- 
coupling took place at densities between 10^^ and 



10^^ cm Further refinements in the 



chemistry 
ram size dis- 
199lh yielded 



(jUmebavashi fc Nakano I |1990D and the 
tribution (|Nishi. Nakano. fc Umebavashi 
similar results. 

Current-driven instabili ties may cause anomalous 
magnetic-flux destruction (|Norman fc Hevvaerts I [r985f ) 
above densities lO^-'^cm"^. Such instabilities, however, 
require drift velocities of charged species greater than 
their thermal velocities, a condition not realized in any 
of th e dynamical calculations ([Fied ler fc Mouschovia£l 
W95, 1993^, 'Ciolek fc MouschoviasJ 1 199,1 119941 119951 : 
Basu fc Mouschovias .1994) . 

Recent dynam i cal calculations by 

iDesch fc Mouschovias I ()2001[ ). which follow the evo- 
lution up to densities n = 2 x 10^^ cni^"^, show that 
ambipolar diffusion is more effective than Ohmic dissi- 
pation in decoupling the magne tic field from the neutral 
matter. [Earlier calculations (jNakano fc Umebavashi I 
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T986allbl: lUmebavashi fc Nakano I [l990l : 

Nishi. Nakano. fc Umebavashi I Il991h overestimated 



the e — H2 elastic collision cross section, and hence the 
Ohmic resistivity, b y a factor of 100; see review by 
iMouschoviasI (|1996f )]. 

In this series of papers, we pursue the self-initiated, 
ambipolar-diffusion-controUed star formation to densi- 
ties ~ 10^^ cm~'^, by which opaque hydrostatic protostel- 
lar cores form. It is in this density regime (after the de- 
tachment of ions and before thermal ionization becomes 
important and the increased degree of ionization reat- 
taches the matter to the field) that the field may become 
completely detached from the matter. A significant con- 
tribution to the resolution of the magnetic flux problem 
of star formation is expected at this stage. Whether am- 
bipolar diffusion or Ohmic dissipation is more effective 
in increasing the mass-to-flux ratio can only be deter- 
mined through such detailed evolutionary calculations, 
not through comparisons of timescales. 

In the present paper we formulate the problem in 
terms of the six-fluid MHD equations (neutral molecules; 
atomic and molecular ions; electrons; negative, positive 
and neutral grains). We do not assume that the mag- 
netic flux is frozen in any fluid (in order to account for 
the possibility that all charges become detached from 
the magnetic field). We derive a new expression for the 
generalized Ohm's law, which lies at the heart of our 
approach to the problem at hand. This new formulation 
(1) includes the effects of inelastic coupling between grain 
species, and (2) explicitly distinguishes between flux ad- 
vection, ambipolar diffusion, and Ohmic dissipation. We 
show how our formulation is related to and can be trans- 
formed into the traditional, "directional" formulation of 
the generalized Ohm's law, and we derive formulae for 
the perpendicular, parallel and Hall conductivities en- 
tering the latter, which include, for the first time, the 
effect of inelastic collisions between grains. In addition, 
we present a general (valid in any geometry) solution for 
the velocities of charged species as functions of the ve- 
locity of the neutrals and of the effective flux velocity 
(which can in turn be calculated from the dynamics of 
the system and Faraday's law). The last two sets of for- 
mulae can be adapted for use in any general MHD code 
to study phenomena beyond star formation in magnetic 
clouds. 

Our results for a fiducial model cloud and a detailed 
parameter study are presented in two companion papers 
(Tassis & Mouschovias 2006b, c; hereafter papers II and 
HI, respectively). As shown in paper HI, at the high den- 
sities treated here memory of the mass-to-flux ratio of the 
parent cloud is lost, and the evolution of the contracting 
core is essentially independent of the initial conditions 
(for values within the parameter space determined by 
observations). This point is of particular importance, 
especially in light of the recent debate concerning the 
formation processes of molecular cloud cores. (For re- 
views and detailed discussion of the arguments on both 
sides of this debate, see, for example, Elmegreen & Scalo 
2004; MacLow & Klessen 2004; Tassis & Mouschovias 
2004; Scalo & Elmegreen 2004; Ballesteros-Paredes & 
Hartmann 2006; and Mouschovias, Tassis & Kunz 2006 
and references therein). The need for a nonideal MHD 
treatment of the high-density regime of core evolution, 
when the magnetic flux problem of star formation has 



to be resolved, and hence the relevance of this work, is 
independent of the adopted theory of core formation. 

2. MODEL CLOUD 

We consider a self-gravitating, magnetic, weakly- 
ionized model molecular cloud consisting of neutral par- 
ticles (H2 with 20% He by number), ions (both molec- 
ular HCO"'' and atomic Na"*", Mg"*"), electrons, neutral, 
singly negatively-charged and positively-charged grains. 
The abundances of charged species are determined from 
the chemical reaction network discussed in Tassis & 
Mouschovias (2005a, App endix A). For a detailed discus- 
sion of the chemistry, see ICiolek fc Mouschoviasi () 19931 . 
119951 ). and iDesch fc Mouschovias I (|200 iD . Cosmic rays 
of energy > 100 MeV, able to penetrate deep inside 
the cloud, drive the ionization for a wide range of 
densities. However, cosmic rays are appreciably at- 
tenuated for column densities > lOOgcm"^. At even 
higher densities, radioactiv e decays become the domi- 
nant source of ionization (Pncuman & Mitchell ' Hill; 
Umebayashi & Nakano 1981; Uniebayashi 1983). The 
ionization rate ^ is the sum of the cosmic-ray ionization 
rate (Ccr) and the rate of ionization due to natural ra,- 
dioactivit y (Crd). The cosmic-ray ioniza tion rate can be 
described (jUmebayashi fc Nakano I[l980l ) by the relation: 



Ccr = CcR,oexp(-cr/96 gem ^), 



(1) 



where a is the column density of the gas, and CcR.o is 
the unattenuated cosmic-ray ionization rate. Values of 
CcR.o are estimated to lie in the range 10 



-17 



10 



-15 



with 5 X 10~^'' s~^ accepted as a canonical value (jSpitzer I 
Il978f ). The ionization rate due to natural radioactivity 
is dominated by the the isotope because of its long 
half-hfe of 1.25 Gyr. In the present calculation s Crd = 
6.9 X 10-23 s-i (|Umebavashi fc Nakano1ll98G[ ). When 
the isotope ^^Al is considered^ Crd = 1.94 x 10"^^ s 



(jConsolmagno fc Jokipii |[l978l ). UV radiation provides 
an additional ionization mechan ism, but it only affects 
the e nvelope of molecular clouds (jCiolek fc Mouschovias I 
[1991 . In this paper, UV radiation is not taken into 
consideration, as we concentrate our attention to the very 
dense, central parts of molecular clouds. 

For simplicity, we consider spherical grains of 
uniform sizej_ since the chemical model used by 



iDesch fc Mouschovias I (|2001h accounted for grains of dif- 
ferent sizes and showed that the effect of the grain size 
distribution on the evolution of the cloud is minimal. In 
the case of ion neutralization through collisions of ions 
(molecular or atomic) with grains, we assume that the 
ions do not get attached to the grains, but that they 
get neutralized and the resulting neutral particle escapes 
into the gas phase. Thus the total abundance of metals 
as well as the HCO abundance remain constant. Grain 
growth and depletion of elements by attachment onto 
grains are not considered here. 

The axisymmetric model cloud is initially in an ex- 
act equilibrium state with self-gravity and external pres- 
sure being balanced by internal magnetic and thermal- 
pressure forces. We follow the ambipolar-diffusion- 
initiated evolution from mean molecular cloud densities 
to a stage beyond the formation of a hydrostatic proto- 
stellar core. Isothermality is an excellent approximation 
for the early stages o f star formation while the densit y 
is below 10^° cm-3 (| Gaustad I IT961 I Havashi I [T96il . 
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The density at which isothermaUty breaks down depends 
on grain opaciti es. With a value o f 0.039 cm^ for 
spherical grains (iDraine fc Lee "1984') , and 0.05 6 cm g 
for fluffy grains (Osscnkopf fc Hcnn ing I [l994f ). the re- 
quired dust column density for the trapping of radiation 
is reached at densities > 10^^ cm~^. However, the opti- 
cal depth is not t he determining factor for the b reakdown 
of isothermality (jMasunaga fc Inutsuka |[l999l ). Instead, 
the rate of compressional heating must be compared to 
the rate of escap e of thermal radiat i on. With these 
considerations. Desch fc Mouschovias I (|2001f ) found that 
deviations from T = 10 K should occur at densities 
> 3 X 10^° cm~'^ but that the temperature is not ex- 
pected to rise rapidly, since a small temperature gradient 
suffices to allow the energy generated by the collapse to 
escape. To avoid the complicated and computationally 
taxing radiative transfer calculation, we switch from an 
isothermal to an adiabatic equation of state at a critical 
density which we treat as a free parameter. 

3. THE SIX-FLUID MHD DESCRIPTION OF MAGNETIC 
STAR FORMATION 

We extend the mu l tifluid formalism described by 
iTassis fc Mouschovias I ()2005D to relax the assumption 
of flux-freezing in the electrons and to account for the 
breakdown of isothermality at high densities. The equa- 
tions governing the behavior of the six-fluid system (neu- 
trals, electrons, ions, negative, positive, and neutral 
grains) are'^: 



^ + V • (p„t;„) = 0, 



(2) 



5(Pgo +Pg- +Ps+) , V7 / , , ^ n 
Ol +V-(pgoi^go+Pg_?^g-+Pg+'i^g+) = 0, 

(3) 

at 47r 



= -e77c(-E + — X B) + 



Vi 



= eni{E+— xB)+F 



(4) 
(5) 

(6) 



= -eng_ {E+^xB)+ Fg_„ + Fg_g„,inci , (7) 
= eng^ {E+'^xB) + Fg^„ + Fg^g„,inci , (8) 

= -Fgon + -fgog-.incl + -fgog+Jnel , (9) 

(10) 

j = e{niVi + Ug^Vg^ - UcVc - ng_Vg_), (11) 
— = -c(V X E), (12) 
V • g - -47rGpn, (13) 

Pn = /(Pn,r). (14) 

^ A set of equations that comprise a logical entity, or a system 
of equations that need to be solved simultaneously, are numbered 
as a block; e.g., Na, Nh, Nc, etc., where N is an integer. 



„ ^ 47r . 

VxB = — J, 

c 



The quantities ps, n-s and Vs denote the mass den- 
sity, number density, and velocity of species s. The sub- 
scripts n, i, e, g_, g+ and go refer to the neutrals, ions, 
electrons, negatively-charged grains, positively-charged 
grains, and neutral grains, respectively. The symbols g, 
E and B denote the gravitational, electric and magnetic 
fields, respectively. The magnetic field satisfies the con- 
dition V ■ B = everywhere at all times. The definitions 
and derivations of F^n (the per unit volume frictional 
force on species s due to elastic collisions with neutrals) 
and F^^s,inci (the per unit volume force on grain fluid 
7 due to the conversion of dust particles of fluid S into 
dust particles of fluid 7) are discussed in detail in §3 of 
ITassis fc Mouschovias I (|2005[ ). 

In the force equations of the electrons, ions, and 
grains, the acceleration terms have been neglected be- 
ca use of the small inertia of these species; it w as shown 
by [Mouschovias. Paleologou. fc Fiedler I (|1985f) that the 
plasma reaches a terminal drift velocity very fast. Simi- 
larly, the thermal-pressure and gravitational forces (per 
unit volume) have been dropped from the force equa- 
tions of all species other than the neutrals because they 
are negligible compared to the electromagnetic and col- 
lisional forces. The inelastic momentum transfer by the 
electron and ion fluids due to attachment onto grains and 
neutralization are negligible compared to the momentum 
transfer due to elastic collisions with neutrals, and they 
have been omitted from the force equations (O and ©. 

The equation of state is either isothermal or adiabatic 
depending on the density. For densities smaller than a 
critical density Popq, the isothermal equation of state is 
used 

Pn-p„C2, (15) 

where C — (fceTiso/M^^H)^^^ is the isothermal speed of 
sound in the neutrals, the Boltzmann constant, ^ the 
mean mass per particle in units of the atomic-hydrogen 
mass toh, and Tiso = 10 K the temperature of the isother- 
mal stages. For densities greater than Popq, adiabaticity 
is assumed and the equation of state becomes 



7-1 

Popq 



(16) 



where 7 is the adiabatic exponent. For simplicity we set 



7 



5/3 for T < 200 K 
7/5 for T > 200 K 



(17) 



to account for the rotational degrees of freedom of H2 
which only become accessible at high temperatures. The 
temperature T is calculated using Equation (|39p . 

For the early stages of star formation and for the 
isothermally infalling gas outside the optically thick in- 
ner region, the temperatures are much smaller than the 
grain sublimation temperature (« 1500 K), and a mass 
continuity equation ^ for the grains (charged -I- neutral) 
is appropriate. 

Altogether, then, we have a system of 13 equations, 
(l2|)- (fT4|) . which contain 17 unknowns (pn, P, E, B, j, g, 

Vn, Ve, Vi, Vg_, Vg_^, Vg^, Po, Pi, Pg„ , Pg_^ , Pg^ ) . To 

close the system, the densities of electrons, ions and 
charged grains (rio, rii, ng^, and ?T.g+) are calculated 
fr om the equilibrium chemical model [see Appendix A 
of ITassis fc Mouschovias 1 ()2005D ]. 
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4. DERIVATION OF A GENERALIZED OHM'S LAW 

The evolution of the magnetic field in a fluid depends 
on the degree of coupling between the field and the fluid. 
If the magnetic field changes in time, an electric field is 
generated in the fluid, as described by Faraday's law of 
induction^ (eq. [H]). This electric field is then capable 
of driving electric currents through the fluid, which in 
turn generate magnetic fields that tend to oppose the 
original change. If the conductivity of the fluid is infinite, 
there are no dissipative losses of the electric currents, 
and the magnetic field is said to be frozen in the fluid. 
The current generated in response to an electric field is 
described by Ohm's law. A generalized Ohm's law can 
be derived from the force equations of the plasma. 

The force equations of all charged species and neutral 
grains, with the elastic and inelastic collision terms writ- 
ten out, are 

= -en^{E + — X B) + —{vn-v^), (18) 



= eni(£; 



Vi 



X B) 



Pi 



= -eng_ [E 



_5z 

c 



X B) 



Ps- 



X B) 



'gn 



' gn 



(19) 



(20) 



(21) 



(22) 



where 



1 



Pe- 



1 



1 



'^gfjo.incl Pgo ''g-iancl 

1 , Ps+ 1 



Pso ''g-g+,inel 
Pe^ 1 



-)-\(23) 
' (24) 



'''goi.incl Pgo Tg+cincl Pgo ''g+g-,incl 

are the relevant timescales for inelastic collisions between 
grains. In particular, t_ is the timescale for a neutral 
grain to be involved in any inelastic reaction involving 
the conversion of negative grains to neutral grains or 
vice versa, and r_(_ is the timescale for a neutral grain 
to participate in any inelastic reaction involving conver- 
sion between positive and neutral grains. 

Dividing each of the first four of the above equations 
by the mass of the corresponding species, multiplying by 
eTsn (where s denotes the species), multiplying the ion 
and positive grain equations by -1, and adding up the 
resulting equations we obtain 



= -E 



i.g_ / gn 



vA X B 



Note that the sequence of events described here does not imply 
causality. An equivalent picture, consistent with Maxwell's equa- 
tions, would be that the magnetic field changes in response to a 
V X E. 



1 / e'^Ua 



+e (^—UcVe + n-iVi — ng_v^_ 



''g+''^g+) 



-en„ 



.(25) 



In equation (|25p . the third term vanishes because of 
charge neutrality, while the fourth term is the current 
density, j. We define the plasma conductivity (jp by 



e "cTon , e n-iTin , e rig 



'gn 



^g+ ''gn 



Too 



TOj 



(26) 



where the individual contributions to cTp are 



2 

e nQTQYi 



^ p,c 



TO„ 



'P,g4- 



g+'''gn 



(27,b) 
(29,d) 



A mean plasma velocity Vp can be defined as 



fp.c 



Vi 



-V, 



fp l^p 
Then equation (j25p is rewritten as 



g- 



"g+ • 



V 



PxB + ^ 

C (Tn 



where 



Jo 



'go 



^go 



'gn 



''"gn 
T+ 



Jo 

0-n 



'gn 



(31) 



(32) 



- v„ 



'gn 



(33) 

The quantity Jq has units of current density and the 
corresponding term in the generalized Ohm's law (eq. 
[32] ) represents the modification of this law because of 
the conversion of neutral grains to charged grains and 
vice versa through inelastic collisions. 

Equation ([5^ is a different form of the generalized 
Ohm's law from the usual expression 



Er, 



??HJ± X es 



(34) 



(jParks 1119911 ). where E^ is the electric field in the frame 
of the neutrals, || and _L denote directions parallel and 
perpendicular to the magnetic field, is the unit vec- 
tor parallel to the magnetic field, and ryu , r]± and ryn are 
the parallel, perpendicular, and Hall resistivities, respec- 
tively. Equation ([5^ can be transformed into equation 
(|34p . with the resistivities appropriately modified to ac- 
count for the effects of the inelastic couplings between 
grains (see Appendix [X)) . However, equation ([5^ offers 
several advantages in explicitly expressing the physical 
significance of each process, distinguishing between dif- 
ferent mechanisms responsible for increasing the mass- 
to-flux ratio (through Faraday's law of induction) and 
simplifying the numerical implementation of the prob- 
lem at hand. 

First, equation ([32|) calculates the electric field in the 
"lab" frame rather than in a frame comoving with the 
neutrals. Although from a mathematical point of view 
the two formulations are equivalent (the electric field in 
the lab frame is simply obtained by subtracting x B/c 
from the electric field in the frame of the neutrals), the 
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interpretation of the results is aided greatly by using the 
lab frame. In this frame, the neutrals are falling toward 
a center of gravity and the magnetic flux follows, at an 
infall speed equal to that of charged species attached to 
the field lines. Flux freezing would be a good approxi- 
mation if the coupling between attached species and the 
neutrals were perfect, in which case the magnetic flux 
would fall in at the same rate as the neutrals. How- 
ever, under the conditions present in the dense cores of 
molecular clouds, this coupling is imperfect, and the flux 
is "left behind" , in a process known as ambipolar diffu- 
sion. After all species except the under abundant elec- 
trons (which interact weakly with the rest of the species) 
have detached, the inward dragging of the magnetic flux 
with the neutral fluid becomes ineffective and its advec- 
tion speed decreases almost to zero, as shown in Paper 
II. 

If this process is viewed in the frame comoving 
with the neutrals, it appears that the flux has to es- 
cape the fluid element by moving outwards. And 
at late times, when ambipol ar diffusion is resurrected 
(|Desch fc Mouschovias I [200ll ). it would be the respon- 
sibility of the extremely few electrons to transport 
this flux outwards - a process that appears improb- 
able and counter-intuitive to some workers, who then 
suggest alternative mechanisms, suc h as grain drifts 
(jNakano. Nishi. fc Umebavashi I [200l ). ^ Such drifts, 
however, cannot possibly be relevant because all charged 
species except the electrons have detached from the held 
lines by the densities of interest (> 10^^ cm~'^), and the 
contribution of detached species to the drift velocity is 
insigniflcant. 

Second, equation ([34|) expresses the components of the 
electric field parallel and perpendicular to the field lines, 
while equation ([5^ splits the electric field in a more phys- 
ical way; when it is substituted in Faraday's law of in- 
duction (eq. [E]), it yields 

^-Vx(t;pXS)-cVx(^ + ^). (35) 
at dp (jp 

The first term on the right-hand side represents fiux ad- 
vection^ (including the effect of collisions, i.e., ambipo- 
lar diffusion), and the other two terms represent flux 
dissipation due to disruption of currents by collisions 
(i.e., resistivity). The last two terms can be combined 
into a single term if one defines an "effective" current 
jj,ff — j In an infinitely conducting medium, fiux is 

simply transported inwards with a velocity Vp, which is 
the conductivity- weighted velocity of the charged species. 

^ ; Nakano. Nishi. fc Umebavashi I lj2002) claim that 
IDesch fc MouschovTasI il200lf ) did not properly account for 
the effect of grains in calculating the Pedersen and Hall conduc- 
tivities, citing equation (28) in IDe sch fc Mouschovias I (120011 ). 
In fact, IDesch fc Mouschovias I pOOl) did not use equation (28) 
to derive their results. They used equations (7)- (13), which do 
account for the effects of all grain species (although not the effects 
of inelastic cou plings between grain species, w hi ch we re also not 
considered by INakano. Nishi. fc Umebavashi I (2002), but are 
explicitly treated here. For the generalized Ohm's law including 
inelastic coupling between grains, see Appendixl^. Equation (28) 
in Dcsch fc Mouschovias (2001) is a simplification presented to 
assist the physical interpretation of the results and is not used in 
any of the numerical simulations. 

® Here we distinguish between flux advection due to the finite ve- 
locities of charged species and any flux redistribution due to Oh mic 
dissipation. The first term on the right hand side of equation I I35I I 
refers only to the former effect, not to Ohmic dissipation. 



Third, if one uses the electron velocity to advect the 
magnetic fiux rather than the ion velocity, then the Hall 
term in the generalized Ohm's law is absorbed in the 
ambipolar diffusion term (see Mouschovias 1987, eqs. 
[72a] and [76]). Here, we use a similar approach, which 
aims at simplifying the generalized Ohm's law by tak- 
ing the advection of the magnetic flux to occur with 
the conductivity-weighted plasma velocity. In this form, 
the components which are responsible for flux dissipation 
and flux advection are clearly separated. Ohmic losses 
are represented by the j^g term, while flux advection is 
represented by the Vp term^. The effect of ambipolar dif- 
fusion is then quantified by the difference between Vp and 
Vn- When all species detach (in which case Vp becomes 
equal to Vn), ambipolar diffusion ceases to operate. 

5. GOVERNING EQUATIONS IN THE THIN-DISK 
APPROXIMATION 

Numerical simulations of the ambipolar-diffusion- 
induced evolution of axisymmetric, isothermal, molecu- 
lar clouds have shown that force balance is rapidly estab- 
lished along field lines and is maintained throughout the 
evolution of the model cloud, even during the dynamic 
phase of contraction perpendicular to the field lines that 
follows the format ion of a thermally and magnetically 
supercritical core (iFiedler fc Mouschovias I \T992. 199j; 
Desch & Mouschov ias 1 120011). F ollowing the approach 
of .Ciolek fc Mouschovias I (|1993( ). we take advantage of 
these results and model the cloud as an axisymmetric 
thin (but not infinitesimally thin) disk, with its axis of 
symmetry aligned with the 2-axis of a cylindrical polar 
coordinate system (r, 0, z) and with force balance along 
the field lines maintained at all times. The sense in which 
the thin-disk approximation is used is that the radial 
variation of any physical quantity is small over a dis- 
tance comparable to the local disk half-thickness. The 
upper and lower surfaces of the disk are at z = +Z{r, t) 
and z = —Z{r,t), while the radius R of the disk satisfies 
the condition Z. The magnetic field is poloidal and 
has the configuration described in Tassis fc Mouschovias 
(2005, eqs. [17a] and [17b]). The model cloud is embed- 
ded in an external medium of constant thermal pressure 

-Poxt- 

For computational simplicity, we reduce the mathe- 
matical dimensionality of the problem by integrating the 
equations analytically over z, assuming z-independence 
of all physical quantities inside the model cloud. This is 
an excellent approximation since thermal-pressure forces 
smooth out density gradients over length scales smaller 
than the critical thermal lengthscale (|Mouschovias I 
Il991l ). which is comparable to Z{r,t). A detailed deriva- 
tion of the thin-disk equations under the assumption of 
isothermality and flux f r eezing in the ion fluid is given in 
ICiolek fc Mouschovias I (|1993f ) . Here we extend the for- 
malism to account for breakdown of isothermality and of 
flux-freezing in the electrons. 

The continuity equations for the neutrals and the 
grains in the disk are obtained through integration of 
equations ([2]) and ([3|) over z, and have the form given in 

An effective flux advection velocity, representing the apparent 
velocity of magnetic field lines when Ohmic losses have also been 
accounted for, can also be defined and is equ al t o {c/B'^)E X B. 
This effective flux velocity is our vi (see eq. |45j ) and should not 
be confused with Vp. 
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Tassis & Mouschovias (2005, eqs. [18a-c]). In the isother- 
mal regime, integration of the z-component of the force 
equation of the neutrals (eq. [3]), with the requirement 
of force balance along field lines, yields 

p„(r, t)C^ = Pext + ^Galir, t), (36) 

where Poxt is the constant external pressure. In the adi- 
abatic regime, the pressure becomes 



P = n„knT 



-Pn 



-pC' 



P-mp 
T 

J- ISO 



T 



(37) 



so the force balance along the z— axis yields in this case 
T 



PnC' 



J- Id 



(38) 



In the isothermal regime, equation (|36p is used to calcu- 
late the density pn (r, t) from the column density dn (r, t) 
and the external pressure. In the adiabatic regime, 
though, equation p8|) must be solved simultaneously 
with the Pn — T relation 



Pn 
Popq 



T 



To 



pq 



(39) 



to obtain the density p-a{r, t) and the temperature T{r, t) 
in the disk. 

Integration over z of the r-component of the neutral 
force (eq. [3]), using the total (thermal-pressure p lus 
Maxwell) stress tensor (jCiolek fc Mouschovias I [19931) to 
combine the thermal-pressure and magnetic-force terms, 
yields 



dt 



dr 



fnffr + -Fmag,r- (40) 



The expression for the total radial magnetic force (per 
unit area), includes the magnetic tension force (the Br 
terms) acting on the surfaces of the disk, and is formally 
identical to Tassis & Mouschovias (2005), equation (ISg). 
The radial component of the thermal-pressure force in 
the isothermal regime is 



r<2 9(Tn 



where 



C. 



off 



TT , [3Pe.t + (7r/2)G'ag] , 
2 " [Pext + (7r/2)Ga2]2 ' 



(41) 



(42) 



and accounts for the contribution of the radial force ex- 
erted by the external pressure P^xt on the upper and 
lower surfaces of the disk, since the surfaces are not hor- 
izontal. In the adiabatic regime, the thermal-pressure 
force becomes (see Appendix [B] for derivation) 

i'V,r - -Ocff 77- - 0„^wO-nTr: 7?^ > (43) 



' T- dr 



T 



where 



Pcxt + ?G(T^ 



(44) 



The r-component of the gravitational field calculated 
from Poisson's equation for a thin disk has the same form 
as Tassis & Mouschovias (2005), equation (18h). 

The assumption that the model cloud is embedded in 
a hot, tenuous, current-free medium implies that the 
magnetic field in the external medium is irrotational 
iy y. B — Q). This allows the use of a scalar magnetic 
potential, from which the radial component of the 
magnetic field at the top surface of the disk can be ob- 
tained (see Tassis & Mouschovias 2005 eq. [18i]). 

In the thin-disk approximation the magnetic field is 
purely poloidal and from Ampere's law (eq. |10| ) the 
current density must be purely toroidal, and thus they 
are orthogonal to each other. Consequently, the £J-field 
has no component parallel to the magnetic field, as re- 
quired by the generalized Ohm's law (eq. [32j). We may 
then define a new quantity, Vf, which has units of velocity 
and can be used to describe the electric field through 



(45) 



Substitution in Faraday's law of induction (eq. [12]) 
yields the evolution equation for the magnetic field 



OB 

— = VxivfX B). 



(46) 



In the plane of the disk (z = 0), equation (jl^ becomes 



dB 



2:,cq 



dt 



1 9(rgz,eqff,r) 

r dr 



= 0. 



(47) 



Note that the velocity Wf defines a frame of reference in 
which the magnetic flux is conserved. 

The velocities in the r and (p directions of all species, 
except the neutrals, and vi can be obtained from the 
solution of the steady-state force equations of the species 
and Ampere's law (eqs. [5]- [9] and [ID]). The derivation 
is given in Appendix |Cj All the radial velocities have the 
form 

1 



(48) 



where Gg is the indirect attachment parameter of species 
s: for 8s ^ 1 the r-velocity of species s is approximately 
equal to that of the field lines, and the species is well 
attached to the magnetic field, while for 0^ ^ 1 the 
species is detached and is falling in with the neutrals 

{Vs « Vn). 

6. THE DIMENSIONLESS PROBLEM AND METHOD OF 
SOLUTION 

We put the equations in dimensionless form by choos- 
ing units natural to the model molecular clouds. The unit 
of velocity is the isothermal sound speed in the neutral 
fluid, C; the units of column density Ucrcf and acceler- 
ation 27rGiTc,rcf are those of the neutral column density 
and gravitational acceleration at the surface on the axis 
of symmetry of the reference state from which the initial 
equilibrium state is calculated. The unit of magnetic field 
is its uniform value in the reference state, B^ci- Although 
not necessary, we define for computational convenience 
the unit of temperature as the temperature of the isother- 
mal lower density parts of the model cloud Tjso — 10 K. 
The temperature above which the rotational degrees of 
freedom of II2 are excited is denoted by Trot ■ 
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We denote the dimensionless form of the cylindrical 
polar coordinates (r, z) by C) and the dimensionless 
time by r, so that the dimensionless thin-disk equations 
become 



an(e,r) =2pn(e,r)Z(e,T), 

9(XgCrn) lc)[Ccr„(Xg-^^g-,{ + Xg+l^g+,« + XgO?^gO,4)] 



9e 



(49) 
(50) 

-0, 
(51) 



4-l[Pcxt+fT2(^,T)], 

4-^/Vo^'q[Pcxt + a2(^,r)]3/5, 



p < Popq ; 



rp3/2 . 

4-5/Vo^'q[Pext + a2(^,r)]5/^T,;,^/^ 

7^3/2 
Popq^rot ^ P I 

(52) 

1, p < Popq ; 

4-2/5[Pcxt+CT2(C,r)]2/Vopq'/^ 

^3/2 . 

Popq <■ P <- Popq^rot i 



dT 



2 

cff 



Popq^rot ^ P I 

(53) 



= Pp,e + fT„g^ + Fn,ag,c ; (54) 



^2 f^'^n 



~^e«TiO-g^ - CncwCTn(C) 



57^ 



(3Pcxt + (tI) 



'-i2 

-'now 



2 



p < Popq ; 

■ I Popq < P J 

(55) 
(56) 

(57) 



mag,^ - 

1 dZ 



1 f / 9ScoQ 
Sc,cq B^,z - Z '^■"'^ 



+ 2Bc.cq Se,z 



dZ_ 



9d0= / <eVn(e')-M(e,^'), 

Jo 




(60) 



POO 

Be AO = - / deeiBcAC) - B,,i]M{^, n m 

Jo 

dB^^cq , 1 9(^-Bc,oq^f.c) 



0. 



(62) 



Equations ([49l -n) and the equations for the velocities of 
all species (Appendix [C| comprise the dimensionless sys- 
tem to be solved numerically. The dependent variables 
are now dimensionless but in this section only, for ease 
in identification, we have retained the same symbols as 
for their dimensional counterparts. 

The free parameters of the problem are: the initial 
central mass-to-flux ratio ^d,cO = (27rG^/^crc,re/)/Sre/ in 
units of its critical value for gravitational collapse in disk 
geometry, l/{2nG^/'^) (|Nakano fc Nakamura 1978); the 
constant external pressure relative to the central grav- 
itational "pressure" along the axis of symmetry of the 
reference state Pext = Pext/(f Gccre/); and the density 
'^opqi above which the equation of state becomes adia- 
batic. 

The reaction rates of the chemical network and the 
value of th e radius of the dust grains ar e given in Ap- 
pendix A of l^ssisASQUSchovias] (IIQQI) . [The thermal 
(Jeans) length does not appear explicitly in the equa- 
tions because of the choice of the units we are using: 
The unit of length C^/27rGcrc.iof is proportional to the 
Jeans length.] 

A control-volume formulation is employed for solv- 
ing numerically the dimensionless equations that gov- 
ern the evolution of the model cloud. The numerical 
method consists of an implicit time integrator; a conser- 
vative upwind advective difference scheme that incorpo- 
ra tes the sec o nd-or der accurate monotonicity algorithm 
of Ivan Leer I (jl979f ): a second-order difference approxi- 
mation of forces inside a mass zone; an integral approx- 
imation of the gravitational and magnetic fields; and an 
adaptive mesh capable of resolving the smallest natu- 
ral lengthscale (i.e., thermal lengthscale Xr.crit)- Using 
convergence tests, we found that 85 zones in the radial 
direction are adequate to follow the evolution of the core 
to the densities of interest with an accuracy better than 
1%. The algebraic equations for the equilibrium abun- 
dances of charged species are solved at each time step. 
A detailed description of the numerical methods is given 
by Morton, Mouschovias, & Ciolek (1994). 

In order to follow shocks that might arise during 
the evolution, we change the originally central-difference 
scheme for the calculation of spatial derivatives to 
a second- order acc urate one-sided difference scheme 
(Paleologo u fc Mou schovias 19 8^. In addition, tensor 
artificial viscosity (Tscharn uter fc: Winkler 1 1 1979( 1 is in- 
troduced. The expression for the artificial viscosity term 
that is add ed to the RHS of the radial for ce equation ([54)l 
IS given m iTassis fc Mouschovias I (|2005l ). 

Our initial and boundary con ditions are the same as 
those discussed in detail in §5 of lTassis fc Mouschovias I 
(l2005h . 

7. DISCUSSION 

We have presented the formulation of the problem of 
the ambipolar-diffusion-controUed formation and evolu- 
tion of magnetically supercritical molecular cloud cores 
inside initially magnetically subcritical parent clouds. 
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and the core collapse to high densities, past the detach- 
ment of ions from the magnetic field and into the opaque 
regime, in which we approximated the treatment of radi- 
ation trapping using an adiabatic equation of state. We 
have used the six-fluid MHD equations to describe the 
model cloud. The magnetic flux was not assumed to be 
frozen in any of the charged species, in order to allow 
for the study of the possible detachment of the lightest 
charged species, the electrons, which would signify the 
complete decoupling of the magnetic field from the mat- 
ter. 

We have reduced the numerical complexity of the prob- 
lem by integrating the equations analytically along field 
lines, and we have used force-balance along the same di- 
rection to calculate the thickness of the formed magnetic 
disk. 

We have derived a new expression for the generalized 
Ohm's law. It includes the effects of inelastic couplings 
between grains, which not only modify the expression 
for the contribution of each charged species to the to- 



tal conductivity, but also account for the contribution of 
neutral grains to the conductivity (because neutral grains 
become charged part of the time). Moreover, this form 
of Ohm's law, when substituted in Faraday's law of in- 
duction, explicitly distinguishes between advection of the 
magnetic flux by charged species and Ohmic dissipation, 
which involves disruption of currents by collisions. 

In two companion papers, we present the solution of 
the problem for a model cloud with fiducial parameters 
(Paper II) and a detailed parameter study (Paper III). 
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APPENDIX 

GENERALIZED OHM'S LAW: EQUIVALENCE OF DIFFERENT FORMULATIONS 

In this Appendix we present the derivation of a generalized Ohm's law (in its usual directional formulation) from the 
six-fluid MHD force equations, accounting explicitly for the inelastic coupling between grain species. The inelastic grain 
coupling not only modifies the perpendicular, parallel, and Hall resistivities, but also introduces a second dissipative 
term associated with the neutral grain fluid. 

In § 2] we presented an alternative formulation of the generalized Ohm's law, which is written in a maximally 
simplified functional form by determining the appropriate, conductivity-weighted plasma velocity, with which the 
magnetic flux is being advected. In this Appendix we discuss how the directional formulation of the generalized Ohm's 
law is related to this maximally simplified formulation. In addition, we present a general solution of the velocities of 
all charged species as a function of the velocity of the neutrals and an effective magnetic-flux velocity. 

General Solution for the Velocities of Charged Species and Neutral Grains 

We adopt a local cartesian coordinate system such that the x— axis is parallel to j x B, the y— axis is parallel to 
(j X B) X B, and the z— axis parallel to B. We can then write 

B = Be, , j = jj_ey + j\\e, , j x B ^ j±Be^ , {j x B) x B = -j^B'^e.y , (Al) 

where e^,. By and are the three orthonormal unit vectors. In addition, we define an effective magnetic- flux velocity 
Vi by 



-Vf X B ~ E±. 



(A2) 



The physical meaning of this velocity is that it represents the velocity of the boundary of a constant-flux surface. The 
effects responsible for the motion of this boundary can be motion of the fluid clement, flux redistribution between fluid 
elements, or flux dissipation due to Ohmic losses. With these definitions, we can write the force equations ©-([HI) in 
component form in the x — y plane as 



= CJcTon(ffa - Vcy) + Vnx - Vex 
= Worcn(Wc2: - Wfi:) + Vny - Vcy 
O^UJiTin{Viy - Vfy) + Vnx - Vix 
= UJiTin{v{x - Vix) + Vny - Viy 

O^UJgT^niVfy - fg-y) + Vnx ^ Wg- 



PgO Tgn 



(VgOx - Vg-x) 



= UJgTgniVg-x -Vfx)+ Vn 



y '^g-y 



— — {VgOy-Vg-y) 

Ps- ^- 

O^UJgTgn{Vg+y - Vfy) + Vnx - Vg+x + ^^ — {VgOx " Vg+x) 

Pg+ T+ 

= CJgTg„(Vf^ - Vg+x) + Vny - fg+y + — — {VgOy " fg+y) 

Pe+ 

T T 

= Vnx - VgQx + — {Vg-x - VgOx) + — {Vg+x " VgQx) 



(A3) 
(A4) 
(A5) 
(A6) 

(A7) 
(A8) 
(A9) 
(AlO) 
(All) 
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= Vn 



VgOy 



gn / \ I ''"gn / 



(A12) 



Once the x— and y— components of the velocity of the neutrals and the magnetic-flux velocity have been specified, 
the X— and y— velocities of the other species can be found by solving the linear, algebraic system of the force equations 
of the species (in steady state). In matrix form, the system of equations (|A3p - (|A12p is 



where A is the matrix of the coefficients of the unknowns 



with 



and 



A= 



To 



" -1 











- 













" 
















-1 

















-1 

























U — 


t^iT in 


U 


U 


n 
U 


U 


1 

1 


U 


n 
U 


U 





- 


-1 - r_ 





r_ 






















'^g 'gn 














— 1 - r_ 





T 











-l-r+ 


r+ 











CJgTgn 














-UJgTgn 














-l-r+ 










''"gn 


'Tgn 























r_ 




To 

































Tgn 


Tgn 


Tgn 
















r_ 




To 




1 

H 




-1 


^ PgO 


Tgn 




_ PgO Tg 








T_ 






Pe- 


r_ ' 


Pg+ T- 


f 





(A13) 



(A14) 



(Al5,b,c) 







" -1 " 









0" 




— l^oTon 


Vix 









^cTon 




-1 







Vg-x 




-1 














t^iTin 


Vg+x 









-t^iTin 




-1 







VgQx 
Vcy 




-1 




, c^^ = 




WgTgn 






-1 




-WgTgn 




Viy 




-1 














WgTgn 


"g-a 

Vg+y 
. '"gOy . 






-1 






-WgTgn 







-1 



-1 










We use Cramer's method to solve the above system. We define 

D = det[A] . 



(A18) 



(A19) 



In addition, we use the notation Z)^^ to represent the determinant of A where the column corresponding to the 
x— component of species s has been substituted by C"^, the notation Z?"^ for the determinant derived from D if the 
y— component of species s column is replaced by C"'^, and similarly for the other C vectors. Then, the solution of the 
system is. 



j~,nx 

SX 

Vsx = ^^V 



D 



jjfx 

^ SX ^, 

nx ^ Vfx 



nny nfy 



for the a;— components of the velocities, and 



jjnx 

sy 

Vsy = -TT-V 



jjtx 

^ ^nx^ V!x 



D 



Jjny 
sy 
-'I 

D 



D 



-Vfy 



^ny 



sy 

D 



Vfy 



(A20) 



(A21) 



for the y— components. 
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The determinants defined above satisfy the relations 

D^}+D%^D6,,, (A22) 

where i and j can take one of the values x, y and Sij is the Kronecker delta. These properties can be easily proved 
as follows. The determinant + D^*- is the same as the determinant which we get if we substitute the sj-column of 

D by + If from this column we then subtract the 3 other j-columns, we get either the determinant of the 

coefficients D {ii i — j), or a determinant which is identically zero (if i ^ j). We define 

= % (A23) 



0s, = -j^ (A24) 

A.. = = -- ^ (A25) 

A.. = -^ = — ^. (A26) 
With these definitions and the properties (|A22[) . the solution can be re- written as 

Vsx = Q Q^"^ ^ ^tx + ^sxivfy " Vny) (A27) 

for the a;— components and 



1 Osy 



T""f + « —-r'"lv + ^sy{vi^ - Vnx) (A28) 

for the y— components. 



^ « 

^ sy T' -L ^ sy 



Drift of Magnetic Flux Relative to the Neutrals as a Function of the Current Density 

Let us define Wg ^ Vg — and Wf — Vf — Vn- Adding up the force equations for all species (including the neutral 
grains) we obtain the plasma force equation, 

B2 1 cr^ 1 

& ^ (uJsTsn) <^P C 

Taking the x— and components of equation (|A29|) we find 

B 



-^pEt r^—Ws = -jxB. (A29) 



"^pEt — —^—'^sx^j±, (A30) 

C ^ [LUsTsn) O"? 

and 

^ (Wsrsn)-^ CTp 

In addition, equations (|A27P and (jA28P give: 

Wsx = 'E^~7'^t^ + ^sxWfy (A32) 
^sx ~r i 

and ^ 

Wsj, = — ^Y^Wfy + AsyWix ■ (A33) 

^ sy ~r 

Hence, equations ()A30P and (|A3ip can be rewritten as 

El CTs Qsx , 1 C /AO.\ 

^ (WsT^n)-^ (Tp Gsx + 1 (WsT^n)-^ (Tp B dp 

and 

^^y E - + "^^y E - - ■ ( A35) 

^ (t^sTsn)^ (Tp " ^ ^ (t^sT,n)2 (Tp G^.y + 1 

For convenience and brevity, we define the coefficients of the above equations to be 

^ {uJsTsny CTp 9si + 1 ' * ^ (WsTsn)^ CTp (A36) 

with i again being either x or y. Then, the solution of equations (|A34[) and (jA35[) with respect to W[x and W[y is 

c ji- (^y c j± Xy (A37) 

B (T-p O^Oy ~ Xx\y B Up OxOy — XxXy 



11 



Equivalence of Equation \32\) and 1^34^ for Ohm 's Law 
Adding (fn/c) x B to both sides of equation ([5^ we obtain 



1 



-Vn X B = - 



— X B 



(A38) 



where is the electric field in a frame comoving with the neutrals. We rewrite the first term in the right-hand side 
of equatio n (IA38|) in terms of Wfx and Wfy, using the definition of Vp (eq. [31j ) to introduce the drift velocities w and 
equations (|A32p and (jA33|l to eliminate Wcx and Wcy in favor of Wfx and Wfy: 



Vp-Vn „ 1 

— X B^ - 

c c 



X B 



= - — (^fc ~ ""n) ^ ^ 

= -} — Wk X B 

C ^ (Tn 



= — 2^ Wky ~ Sy— Wkx 



. B 
-^ex — 
c 



Wfy 



ECTk Qky , CTfc 
, CTn Bi,„ + 1 ^ CTn 



c 



fey 



E 



<Tp Sfex + 1 



(A39) 



Note that k counts over charged species only, while s counts over all species (including neutral grains). Now substituting 
W[x and Wfy using the solution (jA37p . and defining the sums above to be 



(A40) 



where i = x,y, we obtain 



X B^e 



J-L 



_e hi + A ^ 



31- 



e. 



9xOy Aa;Ay 



A, 



9xOy Aa^Ay 



= —J± X 



1 

J 



Qx- 



9xOy A^jAy Qx^y Aa;Ay 



A, 



^xdy ~ '^x\j 



XxKj 



(A41) 



or, defining 



and 



??adv,H 



A 



y 



■'y 



'^X^y A3; Ay OxGy Aa;Ay 

A. 



^x&y — A^jAy 



A.^ 



^y 



(A42) 



(A43) 



we can write the advection term as 



X B = r]Mv,lij± X Bz + ?/adv,±j_ 



(A44) 



The other two components of equation (|A38p can be analyzed as 



CTn CTd CTd 



(A45) 
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and, with the definition Ugo — e^rigoTgn/m, 
Jo _ Joll ,15 1 



(In C CJffTjr 



<J^oB 1 



(Tp C WgTgn 

B 1 



Wfx 



''"gn 



'gn 



'gn 

r+ y e 



^g+- 



+ 



'gn 
Tgn 



e 



g+a; 



gOx 



r- eg_:r + 1 



C'gO 



ZgQx 



Wfy 



^g-\-x ^g—x 



'gn 
T+ 



e 



gOy 



e 



g+y 



®g-i/ 



e 



g+y 



A, 



'gn 



(A46) 



If we now define 



1 (T, 



go 



g 'gn 



9xOy A^jAy 













AgOa; 





'gn 
li 



QgOx + 1 

g+x 



e 



g+a: 



'gn 



e 



T+ Og+x 



-Ag_j; 



(A47) 



and 



1 0-gO 1 



fJp (Tp O^gTgn 

A, 



^x^y Aa;Ay 



Tgn Tgn \ . Tgn . A 
— AgOy + —Ag+y - — Ag_y 

r_ r I / T I T_ 



^X^y Aa;Ay 



'gn 



:0y 



gOj/ 
Tgn 0g+y 



:0a 



e 



we can rewrite the j'g term as 



Jo ^ Jojl 



Substituting equations (|A44|) . (|A45p . and (|A49|) in equation (|A38p . we find that 



^n=^(. 



1 

1- f?adv J 

CTn 



'70, 



'adv,H + '?OHj J_L ^ • 



(A48) 
(A49) 

(A50) 



Note that the inelastic collisions between grains enter this expression both explicitly (through the new term jg|| and 
the new contributions to the perpendicular and Hall conductivities 770^ and ?/oh) as well as implicitly, by altering the 
expressions for r^adv.x and ryadv.H- 

In this Appendix we showed how a generalized Ohm's law can be derived from the six-fluid MHD force equations with 
the inelastic coupling between grains taken into account. We discussed two different formulations of the generalized 
Ohm's law. In the first case, the electric field is separated directionally, into a component parallel to the magnetic 
field, a component parallel to the component of the current which is perpendicular to the magnetic field, and a third 
component perpendicular to both magnetic field and current. In the second case, the contributions to the electric field 
are separated into an advection term and a term quantifying the contribution of each charge species to the Ohmic 
losses when substituted into Faraday's law of induction. The formulations are mathematically equivalent (and we have 
explicitly showed how one transforms into the other). Each offers different advantages for physical understanding and 
interpretation of the formalism and for application in numerical simulations. 

We have presented new formule for all relevant resistivities for the case that inelastic grain couplings are taken 
into account. We have shown that the inelastic grain coupling not only modifies the perpendicular, parallel and Hall 
resistivities, but also introduces an extra dissipative term associated with the resistance to the motion of the neutral 
grain fluid. Finally, we have presented a general solution of the velocities of all charged species as a function of the 
velocity of the neutrals and an effective magnetic-flux velocity. The expressions presented in this paper can readily be 
used in non-ideal MHD codes for the study of weakly ionized astrophysical fluids. Depending on the problem and the 
code, one or the other form of the generalized Ohm's law can be used. 

DERIVATION OF THE THERMAL - PRESSURE FORCE IN THE THIN-DISK APPROXIMATION 

We define a volume element V of the disk as the cylinder-like volume that an arbitrary surface A in the midplane 
(z — 0) sweeps when it is translated vertically until it reaches the top z = Z{x,y) and bottom z = —Z{x,y) surfaces. 
The total force exerted on that volume element V in the i direction is 



F,,^ / V-TdV, 
'v 



(Bl) 
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or, using the divergence theorem 



Fi = j) dS T-n, (B2) 

where T is the thermal-pressure stress tensor, S is the boimding surface of the vohmie V, and n the local unit vector 
normal to the surface. This force can be expressed in terms of the area A on the niidplane, perpendicular to the field 
lines: 



-L 



dAFA. (B3) 

A 

The integrand Fa is the force per unit area on the midplane of the disk. We need to determine Fa- 

We split the bounding surface S into three surfaces: top, bottom, and side. Then the contribution of each one to F 
can be calculated separately, and 

F = Ftop + -^bottom + Fside- (B4) 

The unit normal vectors of the top and bottom surfaces of the disk are 

where V|| is the derivative in the direction parallel to the midplane, perpendicular to the field lines inside the disk. 
The surface element of the top/bottom surface is related to the surface element of the midplane area A as 



dS. 



'top/bottom 



^l + (V||^)2dA (B7) 



The intersection of the side surface with the midplane is a closed curve I. Thus the surface element of the side surface 
is 

dS-side = dl ■ dz, (B8) 
and the normal unit vector is always on the midplane. So 

-Ftop = / dStophtop ■ T{z = Z + e) 



■I' 



dA{z-V\\Z)-T{z = Z + e). (B9) 
Similarly, 

i^bottom = j dA{-z - V||Z) • W{z = -Z- e), (BIO) 

and 

F^ide^ (j)dl J dzi^-W . (Ell) 
The contribution of each surface to the thermal-pressure force is 



-Ftop = 



■ j rfA(5-V||Z)Pext, (B12) 
i=^bottom = - j dA{-Z - V,iZ)Pext, (B13) 



^side ~~ ^ '^^ j dzn\\P, 



z 

dAV\\2ZP. (B14) 



The total thermal-pressure force in the adiabatic regime is 



Fp = -JdA V|| (^(TnC^T^ - 2Pext-^ ) • (B15) 



Force balance in the 2;— direction implies 

p.C2;^=Pext + ^, (B16) 
J- iso ^ 

so that 

(Tn Onl-^ m — CTnLy 7p — 

7 " ^ ISO ^ ISO ^V{^ 7^ 
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Therefore, equations (jBlSp becomes 



T 
T- 

-' ISO 



7rGcr,2 



77 77 P _L £./^/t2 

ISO ISO J ext ~r o n 



rVllCTn 



We define 



and 



_ 2 ' 

"-"off — '-^ 



Then the pressure force is written in a compact form as 



Fp 



dA 



T T 

CnowO'nVll-^ h G^ff— V||(T„ 

-i- ISO -i- iso 



(B18) 
(B19) 

(B20) 



In the r— direction, 



and 



VllCTn = 



dr 



T 



d f T 



ISO \ ISO ^ 

so that the radial thermal-pressure force per unit area in the adiabatic regime is 

F - T da^ „ ^ f ^ 



(B21) 



Clearly, if T = Tjs o the second term vanishes, and equation (|B21[) reduces to the expression (28a) of 
ICiolek fc Mouschovias I (|1993f ) for the thermal-pressure force. If T > Ti^o , the temperature gradient contributes 
to the thermal-pressure force. 



SPECIES VELOCITIES 



Once the velocity of the neutrals, v^, has been specified, the velocities of the other species can be found by solving 
the linear, algebraic system of the force equations of the species (in steady state) and Ampere's law (eqs. [18] - [22] 
and [TUj , where j in eq. [TOj is given by eq. [TT] ) . 

We use equation (|45p to eliminate the electric field from equations (fT5|) - (|2ip : 



0= (Ve -Vf) X B ^ (t)„ - Ve 

C Ton 



{Vi - Vf) X B + —{Vn - Vi) 
C Tin 



en„ 



c 

eng+ 



Vf) X B+^{Vn- Vg-) + — (VgO - Wg- 
Tgn r_ 



C Tn-ri 



ygo 



■"&+)■ 



(CI) 
(C2) 
(C3) 
(C4) 



In addition, we eliminate E (using eq. [45]) and j (using eq. [TO]) from equation p2|) . the generalized Ohm's law, to 
obtain: 



— X B = — 

c c 



1 /cr 



'p,c 



'P,i 



Vi 



cV X B 



T+ 



X B 



(C5) 



where we have also used equations (|5T|) and to eliminate t»p and jo in terms of the velocities of the species. 
Equations (jCl|) - (jCSP and the neutral-grain force equation, 

= — (^ri - -ygo) + — (%-h - -i^go) + — (^g- - ^^go) , 

in component form constitute a 12 x 12 linear system of equations, with unknowns the r— and </>— components of the 
velocities of electrons, ions, neutral, positive and negative grains, and the r— and 0— components of t?/. 



(C6) 
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We now integrate equations (jCip - (jC6|) in the z-direction using the one-zone approximation (no variation of quantities 
with z inside the disk) and decompose them into r— and (/)— components (assuming we have no bulk rotation of the 
disk, Vn<i> = 0). Note that ratios of densities, p, equal ratios of column densities, cr, and the former rather than the 
latter are used whenever such ratios appear. Also, note that the magnetic field contribution in all equations other than 
the generalized Ohm's law (through the gyrofrequencics of the species) only comes from the component. Although 
the integration over z will in principle bring out the contribution of Br from the upper and lower surfaces of the 
disc, such contribution does not actually appear in the ^—components of the equations since it would be multiplied 
by a difference of z— velocities which identically vanishes (there are no drifts in the z— direction). In the case of the 
integrated r— component of the generalized Ohm's law, Br contributes through Fmag- The resulting 12 x 12 linear 
system (with unknowns the velocities fi, Wi^, Wc, i^c^, t'g_ , Dg_0, Wg^ , Ug^^, Ugo, Ugo^, Wf , Wf^) is: 



= UJcTcniVc - Vf) - Vccp 
O = t^iTin(Wi0 -Vf^) + Vn - Vi 
= WiTin(wt ~Vi) - Vi^ 

O = Wgrgn(uf0 - Vg_^) +Vn- Vg_ + 



= 



Pg- T- 
PgO Tgn I s 

Pg- ^- 



O = t^grgn(Wg+0 - Vi^) + W„ - Wg+ + — ^(Wgo - 'fg+) 

PgO 7"gn , 



= CJgTg„(uf -Ug^) -Wg 

= w„ -i^go + ^(wg_ - 



+ 



Pg+ 



(%,0 - Wg+</.) 



£11 / \ ' gn / 

+ (Wg_0 - Wgo</.) + [Vf 



Vi<p = Wc0 H ^10 H Vg_ 4, + 

CTn CTd CTd (T, 



go^ 
+ - Wgo</.) 

p,g+ , 



-^^g+0 



1 



1 



Vi = We + 



cjgr+ 



1 



Wgr+ 



^p 
1 



LJgT+ 



Wgr+ 



(C7) 
(C8) 
(C9) 
(CIO) 

(Cll) 
(C12) 
(C13) 
(C14) 
(C15) 
(C16) 

(C17) 

(C18) 



where tTp.gQ has units of conductivity and is defined as 



'P,go 



(C19) 



and Tpii is a properly averaged coUisional timescale between the plasma (charged species) and the neutrals and is 
defined by 



I ^ Pc (t^cTcn)^ ^ Pi (c^jT-jn)^ ^ Pg+ (^gTgn)^ , Pg-(t^gTgn)^ 



' pn 



Pn 



Pn 



Pn 



Pn 



(C20) 



'gn 



Note that in equation (jClSp . cr„ is the column density of the neutrals, not a conductivity. 

We solve the system of equations (jC7p - (IC17I) in terms of and V[, and the n sub stitute all velocities in equation 
(jClSP to find V{ in terms of Vn- In matrix form, the system of equations (|C7P - (|C17P is 



Vf 



cf 



(C21) 
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where A is the matrix of the coefficients of the unknowns 



A- 



■ -1 














— WeTen 














WeTen 
















-1 




















-1 














WiTin 











-WiTin 





-WiTin 














-1 




















-1 -r- 





r_ 








-WgTgn 








CJgTgn 






















-1 — 'r_ 





T 














-l-r+ 


r+ 











WgTgn 





-t^gTgn 











-WgTgn 














-l-r+ 













Tgn 

r_ 


Tgn 
T+ 


Tgn 
To 









































Tgn 
T_ 


Tgn 
T+ 


Tgn 
To 











V- 


-v+ 


^0 








'^P,g+ 





-1 


0-p 




CTp 


CTp 



and 



To 



PgO Tgn 



r+ = 



PgO Tgn 
Pg+ T+ 



C7DWeT_ 



(7pWgr+ ■ 



(^0 = - </3_ 







- -1 - 




- - 


Vi 









WeTen 


^e- 




-1 
















-WiTin 






-1 


, cf = 












WgTgn 






-1 







fg_,/, 









-WgTgn 






-1 



















- t^f,^ - 




. 








(C22) 

(C23,b,c) 

(C26,e,f) 



(C27) 



We use Cramer's method to solve the above system. We define 

D = det[A] , (C28) 
and use the notation D" to represent the determinant of A in which the column corresponding to s (one of the 
e, e^, i, i^, g_, g_^, g+, go, go0, f^) has been substituted by C". Similarly, D( is the determinant of ^ with the 
column s having been replaced by . Then, the solution of the system is 



dL 



-I 

D 



D 

dL 



-Vf, 



D 



-Vi, 



D 

D 
D 



-Vn + 



D 



-V{, Vg 



D 



g+ 



D 



Vi, V, 



S+1 



D 

D 



-Vn H TT^Vi, 



D 



D 



Di 



gO , -^gO 
■^Vn + -^Vi, Vg„^ 



g+<t> 
D 

j~)n J~)f 



Vi, 



D 



D 

f 



D 



V{,4> = —Vn + -^Vf. 



(C29,b) 
(C31,b) 
(C33,d) 
(C35,f) 
(C37,h) 
(C39) 



Note that for all r— components of the velocities, Dl + D'^ = D. This can be easily seen if in the determinant D( + 
(which is the same as the determinant we get if we substitute the s— column by + C'^) we add to the s— column the 
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3 other columns that also correspond to r— velocities. The result is always the original determinant of the coefficients, 
D. This is a result of the property of the system of equations in which r— velocities, whether known or unknown, 
always appear as velocity differences in the equations. This is not the case for the 0— components of the velocities, 
where the assumption u„0 = breaks the corresponding symmetry. This result is the mathematical expression of the 
physical property of the r— velocities of species, namely, that they vary smoothly between the velocity of the field 
lines and the r— velocity of the neutrals. When a species is well attached to the field lines, its r— velocity is equal to 
that of the magnetic field lines. When collisions with the neutrals become frequent enough for the species to fully 
detach from the field lines, its r— velocity is equal to that of the neutrals. 

We can further emphasize this physical picture by defining ~ {D( + D^)/D, 8^ = /-D", and noting that, as 
we explained above. As = 1 for all r— velocities. We can then rewrite the solution of the system as 



1 



1 



1 



Vf, Vc4, = Ae 



1 



e, + 1 



1 



6,4 



1 



A, 



1 



e 



e0 



e 



e0 



1 



Vf 



Q^4 



e 



i4> 



1 



e 



g+ 



e 



9+ 



Vf, 



■'9+ 

e 



90 



e 



90 



1 



e 



90 



1 



Vf, 



Vfcl, 



A 



/0 



1 



= A 



9+0 



Vf 



a-<t> 



A 



900 



.©9-0 ■ 


f 1 


1 




.Qg+0 - 


— I 

hi 


1 






— V 



e 



-Vf 



9-0 ■ 



e 



9+0 



e 



9+0 



-Vf 



e 



900 



e/0 



1 



e 



900 



1 



Vf 



(C40,b) 
(C42,d) 
(C44,f) 
(C46,h) 
(C48,j) 



Vf 



For the r— velocities Qs is then the attachment parameter (i.e., for 0^ 3> 1, w, 
field lines, while for 0s <^ 1, ~ Wn and species s is detached and comoves with the neutrals). 
Finally, Vf can be obtained as a function of w„ by substituting equations (jCk ) -jCf) in equation (|C18[) 



(C50) 

Vf and species s is attached to the 
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